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Calculation  of  Ray  Paths  in  the  Ionosphere 
Using  an  Analytic  Raytracing  Technique 


1.  INTRODUCTION 


Raytracing  in  a  physically  realistic  model  ionosphere  has  generally  been 
approached  bj  numerically  integrating  the  differential  equations  of  motion  (such  as 
the  Haselgrove  Hamiltonian  equations)  that  describe  a  propagating  ray.  The  method 
developed  here  starts  from  the  paraxial  approximation  to  the  scalar  wave  equation. 
Making  the  paraxial  approximation  allows  the  scalar  wave  equation  to  be  written  in 
the  form  of  a  time  dependent  Schroedinger  equation  where  time  is  identified  as 
distance  along  a  great  circle  path.  Neglecting  diffraction  effects  and  making  the 
ray  approximation  further  simplifies  the  description  to  Newton's  equations  of 
motion.  The  problem  of  ray  tracing  in  a  range  dependent  ionospneric  potential  is 
thus  equivalent  to  solving  Newton's  equations  of  motion  in  a  time  varying  potential 


and  can  be  treated  using  methods  developed  for  a  certain  class  (those  soluble  by 
canonical  transformation)  of  time  varying  potentials.  *  Upon  approximating  the 


variations  in  the  ionospheric  potentials  such  that  they  belong  to  this  class,  analytic 
solutions  for  ray  paths  may  thus  be  obtained. 


(Received  for  publication  18  July  1986) 


1.  Ray,  J.  R.  (1982)  Exact  solutions  to  the  time-dependent  Schroedinger  equation, 
Phys.  Rev.  A,  26:729. 


■Maas 


The  computer  program  described  in  this  report  uses  known  solutions  to  the 
equations  of  motion  for  the  Morse  potential,  which  is  used  to  model  the  E-F9 
ionospheric  potential  well  and  for  the  linear  potentials,  which  are  used  to  model 
the  ground-E  layer  potential  well.  By  fitting  the  complete  range  dependent  ionos¬ 
phere  along  the  great  circle  path  using  shift  and  scaling  functions,  solutions  for 
the  complete  ray  path  are  obtained. 

2.  OUTLINE  OF  THE  PROGRAM 

The  analytic  ray  tracing  program  (RAY1)  consists  of  two  main  parts:  Part  I 
computes  the  parameters  for  fitting  a  Morse  (or  harmonic  oscillator)  potential 
between  the  E  and  F 2  layers  and  a  linear  potential  between  the  ground  and  the  E 
layer  using  previously  generated  ionospheric  data  stored  on  a  file  (TAPE8).  TAPE8 
contains  the  E  and  F ^  layer  values  for  electron  density  peaks,  half  widths,  and 
heights  at  equally  spaced  range  intervals  (usually  1°)  along  a  great  circle  path 
starting  at  the  transmitter  (which  may  be  on  the  ground  or  elevated].  The  results 
of  fitting  the  Morse  and  linear  potentials  in  Part  I  (DO  loop  101  are  stored  in  the 
vectors  RHO(I),  SIG(I),  VMIN(I),  and  SIG(I)  where  I  counts  the  number  of  range 
intervals  between  1-  1  at  the  transmitter  to  1=  MAXR  at  the  receiver. 

In  Part  II,  using  the  series  of  potential  wells  generated  in  Part  I,  trajectories 
are  plotted  for  rays  launched  at  altitude  XSTART  for  a  series  of  takeoff  angles 
(DO  loop  200).  The  ray  trajectories  for  each  initial  angle  TIIET  are  stored  in 
X(J)  and  the  time  delay  is  then  calculated  as  the  quantity  DLAY. 

An  outer  loop  encompassing  loop  10  and  loop  200  (DO  loop  12)  steps  the  fre¬ 
quencies  from  6  to  28  MHz.  After  all  of  the  frequencies  have  been  stepped  through, 
an  ionogram  may  then  be  plotted  by  calling  DLAY  PL. 

To  produce  an  ionogram  for  a  different  path  or  time  of  day  or  season,  a  new 
TAPE8.  DAT  file  must  be  produced  by  running  ION.  COM.  ION.  COM  generates 
the  ionospheric  parameters  along  the  path  selected  using  the  model  ionosphere 
generated  by  IONCAP  with  the  path  parameters  fed  in  as  data  from  the  file 
PVV.  DAT. 


In  this  report,  words  written  in  upper  case  type  refer  to  items  appearing  in  the 
computer  code  and  in  general  are  the  counterparts  of  items  in  normal  mathe¬ 
matical  notation  (which  here  are  usually  lower  case  or  mixed  upper  case  with 
subscripts). 

File  name  notation  conforms  to  the  operating  system  for  the  VAX,  VMS  3.  5. 
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From  the  expression  (1)  for  the  Morse  potential  it  can  be  seen  that  (x,  z.) 

will  have  a  minimum  at  x  =  o(z.)  equal  to  g(z.)  and  an  asymptote  =  g(z.)  +  y  2/p2  at 

l  i  2  i 

x  =  -  oo.  Normally  all  of  the  parameters  a,  y  ,  o( z.),  g(z.)  and  p(z.)  are  fit  at 
range  =  0  (I  =  1).  For  range  steps  with  i  >  1  only  the  minimum  g(z.)  and  its  posi¬ 
tion  a(z.)  and  the  scaling  factor  p(z.)  are  fit.  If  desired,  a  different  range  fitting 
point  may  be  designated  when  initializing  the  program. 

The  analytic  form  for  the  ionospheric  potential  well  is  generated  by  using  the 
l/2(l-x-e'x) 

Chapman  function  fc(x)  =  e  to  represent  the  electron  density  in  both 

the  E  and  F„  layers  and  by  using  a  term  -l/2(l+2x/R  )  to  approximately  account 
for  the  effects  of  curvature  due  to  a  finite  earth  radius  R  .  Its  value  at  any  given 
altitude  x  [in  units  of  100  km]  can  be  found  by  evaluating  the  function  ymm(x)  which 
uses  the  current  values  for  the  ionospheric  parameters  and  the  wave  frequency  v. 
The  complete  expression  for  ymm(x)  in  terms  of  the  layer  half  widths  y  ,  y^, 
heights  xg,  x^.,  and  frequencies  v at  range  point  z^  is  given  by 


u'2  l/2[l-(x-x  )y  -e"(x"xf)/yf] 

ymm(x)  --■='(1+2X/R  )  + - n  e 

l  v*  1/2(1  -  (x-x  )/y  -e  '(x-xe)/ye] 

6 


(2) 


To  facilitate  fitting  the  Morse  potential  parameters,  various  fiducial  points 

are  established  between  the  E  and  F„  layer  peaks  using  the  one  dimensional  mini- 

2  ^ 

miz’ng  subroutine  ZXLSF. 

1.  The  first  point  is  XEEP,  the  position  where  the  ionospheric  potential 
function  YMM(X)  is  a  maximum  near  the  E  layer  peak.  Since  the  position  of  the 
E  layer  peak  is  always  given  by  IONCAP  as  X^  =  1.  1,  XEEP  will  always  be  a  few 
percent  smaller  than  1.  1.  In  order  to  simplify  calculations  of  ray  trajectories  in 
the  linear  ground  -X  potential  well  in  Part  II,  a  point  X1P1  =  1.  05  is  designated 
as  the  peak  of  the  fitted  potential  YPP{X)  near  X£  and  its  value  there  is  defined 

to  be  YMM(XEEP). 

2.  The  second  point  is  S1G(I)  =  SIC-M  which  is  simply  cj(z.),  the  minimum  of 
the  ionospheric  potential  function  YMM(X). 

3.  The  third  point,  XFF,  is  the  maximum  of  the  ionospheric  potential  function 
YMM(X)  near  the  F2  layer  peak. 


2.  ZXLSF  and  ZXPOW  L  are  part  of  the  1MSL  Mathematical  Library,  IMSL  Inc., 
7500  Bellaire  Blvd,  Houston,  Texas. 


4.  The  next  point  to  be  determined  is  XR  which  is  the  point  above  the  minimum 
X  =  SIGM  where  YMM(XR)  =  YMM  (X1P1)  as  shown  in  Figure  1.  The  values  of  the 
potential  at  these  points  are  defined  as  PEFL  =  YMM(XlPl)  and  PEFR  =  YMM(XR) 
and  they  should  be  equal  to  within  the  accuracy  demanded  in  calling  ZXLSF. 

5.  Under  circumstances  where  the  F 2  layer  is  weak  or  for  higher  frequencies, 
it  is  possible  that  the  F 2  layer  peak  will  be  below  that  of  the  E  layer.  In  this  case 

YR  is  set  equal  to  XFF  and  XEE  is*  defined  as  the  point  where  YMM(XEE)  =  YMM(XFF). 

p 

At  this  point  it  is  now  possible  to  compute  the  initial  parameters  y  ,  a,  p( 0), 

a( 0)  and  g(0)  needed  for  fitting  the  Morse  potential  to  YMM(X)  between  X^  and  Xp^* 

For  the  initial  point  I  =  1  we  may  take  p( 0)  =  1  leaving  one  nonlinear  [a]  and  three 

linear  parameters  to  be  evaluated.  As  the  best  least  squares  fit  will  not  necessarily 

have  g( 0)  =  min  |ymm(x)j  we  will  assume  that  g(0)  along  with  y  4  is  one  of  the 

linear  parameters  to  be  evaluated.  Thus  for  some  initial  choice  for  a  and  CT,  the 

2 

conditions  for  a  least  square  fit  for  the  remaining  linear  parameters  y  and  g(0) 
are: 


6  y  x 


A 

L 


y  (x)  -  y  [e 
Jmm  ' 


2loa(x-o)  _  ^2 


"] 


dx  =  0 


(3) 


— >  J  ymm(x)  ( 

Xt 


a(x-a)  .,2  2  i L  .  a{x-a)  ,A  ,  rl  ,  a(x-cr)  ,.5 

-  1)  =  y  J  le  -  1)  dx  +  g  J  (e  -  1) 


dx 


6 


r  r  ,  ,  2  a(x-a)  .2 

ymm(x)  -  y  (e  -  1)  -  g 


dx  =  0 


j  ymm(x)  dx=  y  2  j  (e°^x  ^  -  l)2  +  g  j  dx  . 


(4) 


Defining 


x 

r 

B  =  J  ymm(x)  dx 
Xi 


B 2  =  J  ymm(x)  (eQ'x  ^  -  l)2 

Xl 


(5) 


(6) 


and 


$ 


11 

.vQv 


a  =  I  (e 
n  J 

Xl 


a(x-a)  .  1)2ndx 


(7) 


we  then  have  the  solutions 


y  =  'a0B2  '  a2Bl)/Det 


(8) 


and 


g(0)  =  (a4B1  -  a2B2)/Det 


(9) 


where 


Det  =  (aQa4  -  a2>  . 


The  problem  has  thus  been  reduced  to  a  two-dimensional  minimization  in  a  -  a 


space.  At  this  point  it  might  be  feasible  to  use  a  multi-dimensional  minimization 

2  3 

subroutine.  This  was  tried  using  ZXPOWL  and  FIT.  It  was  found,  however,  that 


even  with  good  starting  values  the  subroutines  tried  would  wander,  producing  non- 
global  minima  in  the  a  -  a  plane. 

To  go  further  we  can  make  use  of  the  relation  satisfied  by  the  Morse  potential 
at  the  point  X1P1. 


2  .  a(XlPl  -  a)  ,.2  . 

y  (e  -  1)  +  g  =  PEFLO 


(10) 


where  the  solution  is  given  by 


a(XlPl  -  a)  =  in  (1  -  [(PEFL0-g)/y  2] 1 12 


(11) 


Thus  for  an  initial  choice  of  a,  with  X1P1  and  PEFLO  given,  a  value  of  a  may  be 

2 

may  be  found  as  a  function  of  y  and  g(0). 

2 

From  the  least  squares  condition  we  have  a  solution  for  y  and  g  in  terms  of 


a  (and  a).  If  we  assume  that  a  solution  exists,  a  consistent  set  of  the  parameters 


3.  FIT  is  part  of  the  DATAPLOT  Language  developed  by  J.  J.  Filliben, 
Center  for  Applied  Mathematics,  National  Bureau  of  Standards, 

U.  S.  Department  of  Commerce. 
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cr,  g,  and  y  may  be  determined  by  iterating  the  two  sets  of  equations,  Eqs.  (8)  and 

2  2 

(11),  using  a  in  Eq.  (11)  to  determine  g  and  y  and  using  g  and  y  in  Eq.  (8)  to 
determine  cr.  The  four  parameter  least  squares  fit  has  thus  been  reduced  to'  a  one¬ 
dimensional  minimization  as  a  function  of  a.  where  for  each  value  of  a  a  consistent 
2 

set  of  y  ,  a,  and  g  is  found  by  iteration. 

These  procedures  are  implemented  in  the  code  by  using  ZXLSF  to  minimize 
the  quantity  SMINIT  (ALF)  which  is  the  sum  of  squares 


jymm  (x.)  - 


a(x.-a) 


-  1)  +  g 


]!2- 


SMINIT  in  turn  calls  GMINIT  to  get  the  self  consistent  set  of  values  for  y  ,  a,  and 
g,  which  are  found  by  iterating  the  two  relations  (8)  and  (11)  IFLAG  times.  The 
final  values  produced  are  stored  as 

SIG(l)  (=cr),  GIM2(=y2), 

VMIN(l)  (=g),  and  ALF  <=e)  . 


To  obtain  p(z.)  and  a(z.)  and  g(z.)  for  fitting  subsequent  ionospheric  profiles 
down  range,  the  least  squares  sum  in  Eq.  (12)  is  minimized  as  a  function  of  p(z.) 

(  =  RHOM  =RHO(I)).  This  is  accomplished  by  calling  ZXLSF  to  minimize 
SUMSQ  (RHOM)  which  determines  cr(z.)  through  the  relation  (12)  using  as  inputs  the 
new  fiduciary  points  (for  the  range  z.),  the  value  of  YMM(XIPI)  (=PEFL)  and  the 
already  determined  ALF  and  GIM2. 

This  procedure  is  performed  MAXR-1  times  within  DO  loop  10,  thus  fitting 
the  E-Fg  potential  along  the  entire  path  length  by  a  varying  Morse  potential  whose 
variation  is  constrained  such  that  exact  z  dependent  solutions  to  the  equations  of 
motion  are  possible. 

For  the  potential  well  between  the  ground  and  tile  E  layer  peak,  we  have  chosen 
to  fit  YMM(X)  by  a  simpler  linear  potential  consisting  of  a  straight  line  segment 
from  X  =  0  to  X  =  X1P l-SIGL(I)  with  a  negative  slope  of  magnitude  GRAV  =  1/R 
with  initial  value  of  -0.  5,  and  a  straight  line  segment  from  X  =  X1P1  -  SIGL(I)  to 
X1P1  with  a  slope  GRVR  =  (PEFL-VGMIN)/ Yg  having  a  final  value  of  EH  above  the 
minimum  at  X  =  X  IP  1 -SIGL(I). 

The  value  of  SIGL(I)  that  will  accomplish  this  is  determined  from  the  value 
of  PEFE(I)  |  =  PEFL- VGMIN]  as 


SIGL(I)  =  PEFE(I)/(GRVR  +  GRAV)  . 


(13) 


From  Figure  2  we  also  have 


EH  =  SIGL  -  GRVR 


DEH  =  SIGL  *  GRAV 


Repeating  this  procedure  at  each  range  interval  I,  we  will  thus  map  out  a  varying 
ground-E  layer  potential  constrained  by  the  requirements  a(x^)=XlPl-SIGL(I), 


p(z.)=l,  that  allow  exact  z  dependent  solutions  to  the  equations  of  motion  between 


x  =  0  and  X1P1  to  exist.  The  subroutine  PLOTWELL  may  be  called  at  every 
IPLWL'th  range  point  to  obtain  a  plot  of  YMMjX)  and  YPP(X)  vs  X  where  YMM(X) 
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Figure  2.  Plot  of  the  Fitted  Linear  Ground-E  Layer  Potential  Well  ypp(x) 
Versus  Height  x 
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is  the  original  ionospheric  potential  well  and  YPP(X)  is  the  potential  well  fitted  by 
a  linear  potential  between  0<X<X  1P1  and  a  Morse  potential  between  X1P1  and  XFF. 


4.  CALCULATION  OF  EXACT  RAY  PATHS 


In  Part  II  of  the  program  the  values  of  cr(z)  and  p(z)  are  used  to  determine  the 
global  ray  path  x(z)  from  the  relation 


x(z)  -  q(z) 
p(z) 


(15) 


where  r(z)  and  p(z)  have  been  previously  determined  for  the  linear  and  Morse 
potential  wells  respectively  and  x'(z')  is  the  solution  of  the  local  equation  of  motion 


=  -V  ,V(x')  . 


(16) 


~  ~  2  ttx1  2 

Here  V(x')  is  given  by  the  Morse  potential  V^(x')  =  y  (e  -  1)  for  rays  in  the 

E-F2  potential  well  and  by  the  linear  potential 


VG(x') 


-  grav  x'<0 
+  grvr  x'- 0 


(17) 


for  rays  in  the  ground  -X£  well. 

The  local  range  parameter  z'  is  given  by 


0  p2  (r) 

9 

Where  z  is  related  to  the  range  £  =  R  /  dS  by 

°  0 


(18) 


(19) 


with  no<C')  the  refractive  index  at  the  minimum  of  the  corresponding  potential  well. 

The  relation  between  the  potential  and  energy  in  the  global  and  local  frames 
is  given  by 


This  yields 


2  -I 

'  _  r  dr  =  j_  r_l_  1  ] 

'■'o  [p0+roT]2  "  ro  Leo  ^  rozi>J  ' 


T-Vti 


■  I  - - - T+S  — — - ■ 

T  =  *i  -  *i  =  0  iPi  +  r!Tl  0  lPo+  rOT]‘ 


=  rL  [Pl  pj  +  rQ  [pQ  pj  * 


rl  =  (P2  ‘  Pi)/Az] 


this  may  be  written  as 


,  l  A22  AZj 

Z2  I  Pi  P2  +  Pq  P 1 


and  in  general 


i  i  A  z.  , 

'  l+l 

z.  ,  =  z.  + -  . 

l+l  i  p.  .  p. 

Mi+  1 


This  computation  appears  in  the  program  before  label  919  if  the  ray  is  propa¬ 
gating  in  the  upper  well,  for  which  the  indicator  IFLAG  is  set  =  0.  If  the  ray  if 
propagating  in  the  lower  well,  IFLAG  =  1  and  z ,+  ^  is  calculated  at  the  line 


labeled  109. 


For  a  ray  launched  in  the  E-Fg  duct,  the  initial  energy  with  respect  to  the 


bottom  of  the  potential  well  (VMINO)  is  given  as  E  =  -VMIN0  +  VSTRT  +  COS(t) ) 
where  VSTRT  is  the  potential  at  the  point  XSTRT.  Depending  on  whether  E  is 
greater  or  less  than  Uq,  z'  =  ZP(1)  is  set  such  that  the  ray  will  have  positive  or 
negative  velocity  according  to  whether  6  is  greater  or  iess  than  zero.  The  absolute 
value  of  the  velocity  is  that  which  a  particle  released  from  the  right-hand  turning 
point  where  the  line  E  =  const,  intersects  YMM(X),  would  attain  on  arrival  at 
XSTART.  The  initial  value  for  ZP  is  calculated  above  label  1255  an-'  then  used  to 
calculate  the  global  position  X(J)  and  local  velocity  XPDOT  inside  the  range 
DO  LOOP  on  label  20.  These  quantities  are  needed  at  each  range  step  to  evaluate 
the  time  delay  for  the  entire  path  given  by 


ft 


t*. 


m 


w 


where 


d?  =  (Ro  +  x)  de  ,  §  =  10  +  p  X'  +  g;  £  ] 


dx' 


P 


and 


x(z)  =  p(z)  •  x'  (z'(z)]  +  a(z)  .  (38) 

This  is  accomplished  by  using  Simpson's  N  point  integration  with  the  integrand 
calculated  at  label  101  if  X(J)  is  in  the  upper  well  and  at  label  201  if  X(J)  is  in  the 
lower  well.  The  final  sum  is  converted  into  milliseconds  below  label  20. 

For  ray  paths  in  the  upper  well,  situations  may  occur  for  high  angle  launches 
and  at  high  frequencies  where  the  fitted  potential  well  prohibits  a  ray  from  escaping 
over  the  Fg  peak  even  though  it  is  energetically  possible.  To  remedy  this  situation 
a  further  fiducial  point  XFFI(I)  is  calculated  at  label  1001  to  denote  the  point  at 
which  the  line  V(X)  =  PEAKF  +  DVMIN  intersects  YPP(X)  where  DVMIN  is  the 
difference  between  the  minima  of  YPP  and  YMM.  If  a  ray  goes  beyond  this  point 
then  it  is  considered  to  have  escaped,  even  though  it  would  have  been  turned  around 
by  the  fitted  potential.  To  check  whether  this  occurs  it  would  be  simplest  to  check 
each  point  in  the  path  to  see  whether  or  not  it  exceeds  XFFI(J).  This  will  not  work 
in  every  case  as  a  ray  may  be  situated  at  X(J)  just  before  the  turning  point  at 
range  Z(J)  and  be  situated  at  an  essentially  similar  point  X(J+1)  after  reflection 
at  range  Z(J+1).  Both  X(J)  and  X(J+1)  may  be  less  than  XFFI(J)  whereas  the 
turning  point  XTURN  could  exceed  XFFI(J).  To  account  for  this  situation  the  quan¬ 
tity  WSW  =  XPDOT(  J)*XPDOT(  J+ 1)  is  calculated.  If  WSW  is  negative,  it  indicates 
that  the  ray  has  been  reflected  between  steps  J  and  J+l.  If  this  is  true,  XTURN  is 
calculated  and  checked  against  XFFI(J)  to  see  if  the  ray  escaped.  These  operations 
occur  in  the  program  just  before  label  3030. 

A  ray  with  negative  velocity  that  propagates  to  the  left  of  X1P1  without  being 
reflected  is  assumed  to  be  propagating  in  the  lower  well  as  though  it  were  launched 
at  X1P1  with  velocity 


where  z^  is  the  range  point  at  which  the  ray  crossed  X  =  X1P1.  Due  to  the  nature 
of  the  solutions  of  the  equations  of  motion  the  energy  E1  measured  with  respect  to 


the  potential  minimum  g^(z  J  (=VMIN(J))  will  be  an  invariant  as  long  as  the  particle 


remains  in  the  upper  well  and  similarly  will  remain  an  invariant  as  long  as  the 
particle  remains  in  the  lower  well  (including  ground  reflections),  e'  may  change 
if  a  ray  returns  to  the  upper  well  at  range  z  >  z^  after  having  propagated  in  the 
lower  well. 

A  similar  argument  can  be  made  for  E^  being  an  invariant  and  changing  only 
for  a  ray  returning  from  the  upper  well.  These  changes  are  calculated  at  the  entry 
point  into  the  lower  well  (label  3031)  and  at  the  entry  point  into  the  upper  well 
(label  301).  For  rays  crossing  X1P1  in  either  direction,  the  excess  in  z',  that  is, 
that  which  takes  the  ray  beyond  X1P1,  is  calculated  and  converted  to  the  starting 
z'  value  for  the  ray  in  the  adjacent  well.  This  is  done  at  label  1150  for  rays  coming 
from  the  right  of  X1P1  and  above  label  225  for  rays  coming  from  the  left  of  X1P1. 


For  the  case  of  a  particle  coming  from  the  right  of  X1P1  with  E  >  Uq  the 


general  solution  written  as 

ax'  [1  -  (y  /^)21 

e  =  u  =  - - - 1 - 1 - J — 

(y/JW')  [ch(a>b(E')z')- y  He<] 


may  be  solved  for  z'  yielding 


*Je'  y  u 


z  ln  \  -  - 

b  '  [E'  -  y2(l  -  u)+  n/e  '  -  y  2  VR] 


where 


R  =  E'  -  y2(l  -  2u  +  u2)  . 


The  crossing  excess  in  z'  in  the  lower  well  will  thus  be  given  by 
Az1  =  [z'(x'  )  -  z'((X1P1-  ad/p.)]  »p.  •  p._ j  . 


The  general  solution  for  a  ray  in  the  lower  well  is  given  by 


x'  =  z'  (v  +  -r  z'  grav)  x'  s  0 


x'  =  z'  (v  —  z1  grvr)  x‘  >  0  . 
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We  assume  that  the  starting  value  of  z  for  a  ray  in  the  lower  well  is  zero  when 
the  ray  starts  from  x'  =  0  moving  to  the  right  or  left  with  velocity  v  =  ±  vq  =  V2Eq  . 

If  no  boundaries  were  interposed  at  x  =  0  or  removed  at  x  =  X1P1  the  range 
interval  z'  required  to  reach  the  left  and  right  turning  points  would  be 
TG  =  |  VO  |  ,/GRAV  and  TGR  =  |  VO  |  /GRV R  respectively.  The  interval  in  z'  to  go 
from  x  =  0  to  the  left-hand  turning  point  is  given  by 

DTL  =  [2*(+ VGMIN  +  0.  5  +  EG  +  DEH)] 1  ^  2/GRAV  (44) 

and  for  x'  =  X1P1  to  the  right-hand  turning  point  the  interval  is  given  by 

DTR  =  [2(EG  -  EH)]1/2/GRVR  .  (45) 

For  EG  greater  than  (VGMIN  +  0.  5  +  DEH),  the  ray  will  be  reflected  from 

x  =  0  in  classical  ground  hop  fashion.  In  the  program  this  is  achieved  by  adding 

2  x  DTL  to  z'  for  a  ray  with  negative  velocity  whose  path  goes  below  x  =  0. 

For  rays  trapped  in  the  lower  potential  well,  the  motion  is  divided  into  four 

dx* 

quadrants  (denoted  by  IQUAD  in  the  program:  (l)SIGMA  £  x  <  X1P1,  <  0; 

H  x 1  rlx1 

(II)  0  <  x  <  SIGMA,  <  0;  (III)  0  <  x  <  SIGMA,  >  o:  (IV)  SIGMA  <  x  <  X1P1, 
dx' 

gfr-  >  0  where  SIGMA  =  X1P1  -  SIGL(J). 

For  motion  in  the  lower  well  the  calculation  is  directed  at  label  1050  to  the 
quadrant  in  which  it  had  previously  been  propagating  (at  range  z.  j)  and  is  checked 
to  determine  whether  X(J)  is  beyond  the  next  quadrant  boundary;  if  not  it  is  allowed 
to  propagate  according  to 

dx'2  _  _  |  grvr  in  (I,  IV) 
dz1  I  -grav  in  (II,  III) 

If  the  ray  crosses  the  I-II  or  III-IV  boundary  at  x  =  SIGMA  the  excess  range  interval 
DZP  is  calculated  and  used  as  the  initial  range  value  in  the  next  quadrant.  A  ray 
in  quadrant  IV  with  EG  -  EH  will  not  be  turned  around  by  the  potential  before 
reaching  X  =  X1P1  and  will  cross  into  the  E-Fg  potential  well.  The  calculation  is 
then  directed  to  label  301  when  the  ray  can  be  considered  as  propagating  in  the  upper 
(Morse)  potential  well  with  energy 

E  -  (EG  +  VGMIN  +  DEH- VMIN(J)  (47) 

and  with  an  invariant  energy  given  by  Eq.  (20)  or  equivalently  by 


p  • 


(48) 


E1 


dx'  (  a  +  p  x'  \ 

df  V  P  ) 


Once  launched  a  ray  will  propagate  in  one  or  the  other  of  the  two  potential 
wells,  crossing  the  boundary  at  X1P1  when  energetically  allowed.  When 
J  =  MAXR  the  ray  path  is  checked  to  determine  whether  the  ray  could  have  been 
intercepted  by  a  receiver.  For  a  receiver  on  the  ground  this  is  done  at  label  880 
by  checking  whether  the  indicator  for  a  ground  reflection  [VREV(MAXR)  or 
VREV(MAXR-l)]  set  at  label  222  is  positive.  If  positive  the  point  of  reflection 
ZGRND  is  calculated  and  checked  to  see  whether  it  is  within  the  receiver's  antenna 
beam  (taken  here  to  be  100  km).  For  rays  which  could  have  been  intercepted,  the 
delay  time  DLAY  is  stored  in  the  array  ADLAY  (ITHET,  1FREQ)  which  it.  subse¬ 
quently  plotted  in  the  subroutine  DLAYPL  as  soon  as  all  frequencies  in  the  DO  loop 
12  have  been  propagated.  For  a  satellite  receiver  in  the  E-Fg  duct,  an  altitude 
bin  is  set  at  label  870  and  all  rays  whose  final  height  XX(MAXR)  is  within  10  km 
of  the  satellite  height  SXAT  are  accepted  and  stored  in  ADLAY  for  plotting. 
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